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Using exact diagonalization methods we study possible phases in a one-dimensional model of two 
differently populated fermionic species in a periodically driven optical lattice. The shaking amplitude 
and frequency are chosen to resonantly drive s — p transition while minimizing the standard intra¬ 
band tunnelings. We verify numerically the presence of an emergent density wave configuration 
of composites for appropriate filling fraction and minimized intra-band tunnelings. The majority 
fermions moving in such a lattice mimic the celebrated Rice-Mele model. Far away from that region, 
the structure changes to a clustered phase, with the intermediate phase abundantly populated by 
defects of the density wave. These defects lead to localized modes carrying fractional particle charge. 

The results obtained are compared with earlier approximate predictions. 


I. INTRODUCTION 

Ultracold atoms trapped in optical lattices provide sys¬ 
tems characterized by an unprecedented control over var¬ 
ious parameters, enabling a simulation of a wide array of 
exotic solid state models. One example of such phenom¬ 
ena are topological insulators [1, 2] that are of particular 
interest in the field of quantum information and spintron- 
ics due to their inherent stability and transport proper¬ 
ties [3-5]. Lattices hosting systems showing topological 
properties have been realized experimentally, both for 
two-dimensional (2D) [6-11] and one-dimensional (ID) 
models (e.g. SSH/Rice-Mele dimer [12, 13] in [14] or 
Thouless pump in [15, 16]). Optical lattices do not al¬ 
low by themselves for generation of impurities on which 
boundary localized modes may appear - the lattices are 
necessarily perfect. In 2D, the localized defect - a vortex 
- may be created by a vortex wave [17, 18] leading to a 
well placed dislocation. In ID, the situation is not so sim¬ 
ple, but a recent proposition [19] suggests that topologi¬ 
cally nontrivial states may emerge in systems consisting 
of two subspecies of strongly attracting fermions. There, 
the topological structure is not encoded in the underly¬ 
ing lattice geometry, but rather is an emergent feature 
arising from atomic interactions, enabling creation of de¬ 
fects with less constraints. For high enough values of the 
interaction strength, fermions of different species tend to 
bind together forming composites, and if there is some 
imbalance in a number of atoms of both species, excess 
fermions stay unbound. To extract essential properties 
of the system, one has to take into account higher bands 
{p band at least, as in the model studied below) and the 
effects of strong interactions, such as the density induced 
tunnelings [20-23]. The lattice shaking is employed with 
the shaking frequency such that the inter-band density- 
induced s to p tunneling is resonantly enhanced. As a 
result, in a ID chain, the emergent system is proposed 
to be described by the Rice-Mele model [20]. For a tri¬ 
angular lattice geometry, similar processes lead to the 
creation of synthetic non-Abelian fields in an emergent 
dice lattice [24]. 

Let us note, parenthetically, that physics of p-orbital 



FIG. 1. Visualization of the kinetic tunneling processes 
present in the system. Blue (open) and red (filled) half¬ 
circles denote the minority, j"-fermions and the majority, 4" 
fermions, respectively. Observe that the minority fermions 
appear paired in composites only. 


fermions is very rich, leading to a possible creation of 
FFLO states [25] as well as density stripes at appropri¬ 
ate fillings due to nested Fermi surfaces [26, 27] even in 
the absence of any periodic driving (for a review of these 
effects see [28]). Those systems were studied using both 
two- and three-dimensional models. Here, we shall re¬ 
strict ourselves to small ID systems amenable to exact 
diagonalization. 

Let us stress that the main approximation used in pe¬ 
riodically driven models discussed in [20, 24] is to neglect 
the tunneling of the minority components. As a result, 
one generates a modified Falicov-Kimball like model with 
immobile composites (made out of strongly coupled pair 
of fermions) and mobile excess fermions. We test this 
assumption in the present paper. Namely, we are em¬ 
ploying an exact diagonalization method to the system 
described in [20] in order to assess the validity of the re¬ 
sults presented there. The fidelity and structure factor 
analysis allow us to classify the ground states for different 
values of parameters. We consider also explicitly possible 
configurations with a given number of defects. 
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large energy. As tested by us with Wannier functions for 
different lattice depths, |[/i |, the energy of the interaction 
between a fermion in the p and in the s band, is smaller 
than \Uo\. 


FIG. 2. Visualization of the direction dependent inter-band 
density dependent tunnelings in the dressed picture. The 
dressed tunneling amplitudes are direction dependent. 


II. THE SYSTEM 


The system considered is a mixture of two species 
of unequally populated, strongly attractively interacting 
fermions in a ID periodically shaken optical lattice. The 
Hamiltonian of the system is [20] H = ^tun + ^dit + 
Hons + Hsh{t), where; 
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With that assumption, the lowest energy manifold is 
filled with composites - pairs of t and ^ fermions - and 
the remaining ^ fermions, leading to nontrivial dynamics. 
Note that, for example, if a minority ^-fermion tunnels 
from a given site, it leads to breaking of the compos¬ 
ite. It costs a huge amount of energy (|t/o|) unless the 
tunneling occurs to a site in which a majority ^-fermion 
waits to form a composite with the t-particle. Only the 
latter process remains in the low energy manifold. In ef¬ 
fect, the simple tunneling of the minority fermion may be 
viewed as a tunneling of the composite, accompanied by a 
reverse direction tunneling of the majority fermion (com¬ 
pare Fig. 2) within this manifold. The system may be de¬ 
scribed by operators describing excess majority fermions 
(residing either in s or in p band) and the composites de¬ 
scribed by annihilation (creation) operators G(ct) obey¬ 
ing hard-core boson commutation relations. The corre¬ 
sponding composite number operator is = cjc^. The 
presented intuitive picture is fully recovered on a more 
formal level by an appropriate construction of the ef¬ 
fective Hamiltonian [19]. In effect, the excess majority 
fermions move in the emergent lattice created by the 
composites. To avoid excessive repetitions we refer the 
reader to [19] for details while [24] provides yet another 
example of a two-dimensional construction based on the 
idea described above. 


Above and in the following, st, are creation and 

annihilation operators of ^-fermions in the s- and p-bands 
respectively, while are s-band creation and annihi¬ 
lation operators for ^-fermions. Accordingly, and 

n\ are the corresponding number operators. Note that 
while we take into account s and p bands for ^-fermions 
we consider only s band for ^-component. That is so 
because we assume that ^-fermions form a minority com¬ 
ponent with filling close to 1/2. On the other hand, we 
assume a bigger density for fermions. 

The single-particle tunneling part of the Hamiltonian is 
given by ^tun- We assume both species to have the same 
mass and feel the same optical lattice for simplicity. With 
the adopted sign convention Jo, Ji > 0. The density de¬ 
pendent tunneling part is denoted as ^dit- The tunneling 
coefficients To,Ti,Toi are given by appropriate integrals 
of Wannier functions [20, 23]. Since the p-Wannier orbital 
is antisymmetric, the inter-orbital tunneling amplitudes 
have opposite signs in opposite directions as reflected by 
{j — i) factor. 

The basic assumption of the model is that attraction 
between different species dominates the problem energet¬ 
ically. Consider the onsite energy term Hons- Under our 
assumption, Uq is negative with \Uo\ giving the large en¬ 
ergy scale. El - the energy of the p band - is another 


The second important step is to derive the effective 
Hamiltonian valid for the high frequency driving obeying 
the (almost) resonant condition 


Ui+Ei=Nuj^A ( 2 ) 


with N being an integer and a small detuning |A| uj. 
Observe that the time-dependent part of the Hamilto¬ 
nian, ^sh(0, contains two time-periodic terms. The first 
one describes a standard horizontal lattice shaking (after 
an appropriate gauge transformation) as originally pro¬ 
posed in [29] and reviewed e.g. in [30]. Such a horizon¬ 
tal shaking has been realized experimentally by several 
groups [31-33] and serves as a convenient knob on lat¬ 
tice system properties. The second term is due to the 
harmonic variation of the lattice depth. This translates 
into a periodic modulation of the p-band energy offset 
6El [20] . The phase p between the two harmonic modu¬ 
lations can be easily controlled in experiments. The pro¬ 
cedure of averaging is fairly standard and is described in 
detail in [20]. We quote here the final effective Hamilto¬ 
nian expressed in terms of composite and excess fermion 
operators: 
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where the tilde sign over tunnelings and density de¬ 
pendent tunnelings indicates their dressed character (af¬ 
ter time averaging). Explicitly, Ji = Jo (^) Jz with 
I = 0,1 corresponding to s and p band respectively 
and Jo being the Bessel function [20, 29]. A simi¬ 
lar dressing takes place for intra-band density depen¬ 
dent tunnelings T/. On the other hand, the inter-band 
s — p tunneling amplitude value becomes direction de¬ 
pendent due to the phase difference between shaking am¬ 
plitudes. We express that asymmetry by denoting the 
tunnelings between pi ^ as Tq^ and pi ^ as 
ToJ These tunneling processes are visualized in Fig. 2 
and read explicitly [20] Tq^ = Jat (A^/cj) Tqi, where 
= y /{K ± 6El cos ^ K‘^ sin^ p. 

While the frequency of the periodic drive is fixed by 
the resonance condition (2), the shaking amplitude K 
provides a convenient parameter to tune the properties 
of the system. In particular, Kc such that KdoJ ~ 2.405 
corresponds to the zero of Jo Bessel function. For such 
a choice of K, the intra-band tunnelings almost vanish 
and the inter-band density dependent tunneling becomes 
the only mechanism of transferring the majority fermions 
(the composites becoming immobile in this limit). Then, 
as suggested in [20] for = 1/2 (the filling for minority 
fermions), the composites form a density wave (DW) in 
the ground state while the excess majority fermions are 
described by a Rice-Mele topological dimer model. On 
the other hand, sufficiently far from Kc the standard tun¬ 
neling mechanisms dominate - the system then organizes 
into a clustered phase (CL) with composites and empty 
sites separated in space [20]. 

To test this prediction, one has to carefully estimate 
various parameters appearing in the minimal Hamilto¬ 
nian, (3). They depend on the details of the lattice po¬ 
tential and interactions between two species. We fol¬ 
low the assumptions of [20] and assume the optical lat¬ 
tice potential to take the form Fiatt = V\\ sin^(7rx/a) + 
FL(sin^(7r^/a) + sin^(7r2;/a)), with a being the lattice 
constant. For V± ^ the system is effectively one¬ 
dimensional. We take V± = 25 while hji = 8 in the units 



FIG. 3. Ji, Ti and tunnelings for Vb = 8, W = 25, 
a — —0.1, iV = 1 and 0 = 0. For K — Kjoj range shown 
in the figure there exist three phases: the density wave (DW) 
(blue), the clustered phase, CL (white) and the mixed phase, 
MX (pink). Boundaries of these phases were calculated using 
the fidelity susceptibility (see fig 4) and the structure factor 
(fig. 5). 


of the recoil energy Eji = iJliJMcJ) (note that a = Xf 2 
with A being the wavelength of the laser beams forming 
a standing wave pattern). As a dimensionless interac¬ 
tion strength we take a plausible value a = asja = —0.1 
(with as being the (negative) scattering length). That, 
together with lattice parameters, allows us to estimate 
all the tunneling and interaction parameters of the model 
using the Wannier functions appropriate for the lattice 
potential [20]. 

As far as the shaking is concerned, we obviously con¬ 
centrate on the vicinity of K/uu = 2.4 region, taking 
the vertical shaking to be in phase with the lateral one 
{ip = 0), which gives Toi > Tq^. For simplicity, we as¬ 
sume first the exact driving resonance A = 0. In Fig. 3 
we show the dependencies of the different dressed tun¬ 
nelings as a function of Kjuj (we shall later assume a 
notation K = K/uj) coming from Wannier function cal¬ 
culations. 

To find the ground state of (3) we have yet to define 
the density of majority component which is taken to be 
unity (thus we have a 1/2-filhng of composites and 1/2 
filling of excess fermions). Then, we use the exact diag- 
onalization method based on [34, 35]. Diagonalizations 
take place in Fock space of all possible configurations of 
the system, assuming that each site i may be empty or 
occupied by a single composite or a 0 fermion in s state, 
or both the composite and 0 fermion, although the second 
one in p state (because there is already s-state 0 fermion 
in a composite). Therefore, the local Hilbert space con¬ 
sists of 4 states per site with no truncation. For an even 
number of fermions, a fermion tunneling between arbi¬ 
trary “edges” (that is, between the first and the last site) 
leads to an additional phase (sign) change arising from 
the anti-commutation relations. Because the number of 
fermions is half the number of sites, available numbers of 
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sites are of the form of L = 4/ + 2, / G Z. 

With periodic boundary conditions, Hamiltonian in 
Eq.(3) commutes with the translation operator (T), 
which allows us to use states with the conserved to¬ 
tal momentum (k) as our basis: T\a{k)) = e'^^\a{k)). 
States with different ks are orthogonal to each other, and 
k G {(—1//2 -f 1)^,..., (I//2)^} (because = 

\a{k)) = \a{k))) with L being the length of the chain. 
Diagonalization consists of creating states in the basis 
(for some/all values of /c), calculating matrix elements 
of H in that basis and using numerical algorithm to get 
eigenvalues and eigenvectors for the lowest energy states. 
We would like to point out that the total momentum 
k serves only to split the large Hamiltonian matrix into 
smaller blocks. 


III. RESULTS 

We carry out exact diagonalizations typically on a 
chain of length L = 14 (leading to matrices of the rank 
^ 840000). For selected data we show the results for 
L = 18 (matrices of rank around 131 million). As tun¬ 
nelings are nearly symmetric with respect to Kc ~ 2.405 
(only T(^ are noticeably different, which leads to small, 
quantitative - but no qualitative - changes), we will 
only consider K < 2.405. In the interval of interest, 
the ground state corresponds to k = n. To character¬ 
ize its properties and locate possible phase transitions 
we use the fidelity approach [36]. We calculate the fi¬ 
delity, F, associated with a small parameter change 6 
(here F{K,6) = F S/2))) using the 

eigenvectors coming from the diagonalization. For ^ 0, 

we get F ^ 1 — which defines the fidelity suscep¬ 
tibility, X [37, 38]. It is commonly understood that the 
fidelity susceptibility diverges at phase transitions. For 
our finite system, the possible crossovers will be identified 
by the maxima of x- 


A. Resonant case 



FIG. 4. Fidelity susceptibility, x? fo^ the ground state for 
L = 14. Observe a sharp peak at A = 2.3 and a rich structure 
of peaks up to about K = 2.35. The different phases are 
analyzed using the structure factor - see the text and fig. 5 



FIG. 5. Structure factor S{K,q) calculated for L = 14 sites 
and different q values as indicated in the figure. For K > 2.35 
a single q — 1 value dominates indicating DW phase. Sim¬ 
ilarly For K < 2.29 q = 1/7 dominates pointing out to the 
clustered phase (GL). Ghanges in the structure factor behav¬ 
ior nicely correlate with peaks in the fidelity susceptibility - 
compare fig. 4). 


As mentioned above, the simplified analysis [20] pre¬ 
dicts the existence of two composite arrangements: the 
density wave (DW) close to Kc = 2.405, where intra¬ 
band tunnelings are effectively switched off, and the clus¬ 
tered phase (CL), where the composites group together. 
Thus, we should expect a single x maximum for K < Kc 
corresponding to the border between these two phases. 
The numerical results are, however, quite different (com¬ 
pare Fig. 4). There are indeed two regions of low fi¬ 
delity susceptibility for K < 2.29 (with a sharp peak of 
X around K = 2.3) as well as for K values close to Kc (for 
K > 2.35) indicating stable phases. On the other hand, 
the interval K G (2.29,2.35) shows a structure of peaks 
with X having significant values almost everywhere. 

To understand that somewhat complicated behavior of 


X we consider the structure factor here defined as 



where fi^ is the number of composite bosons (0 or 1 in 
our case) occupying i-th site. Fig. 5 shows three areas 
with different behavior of structure factor, each corre¬ 
sponding to different phase structure. For the CL phase, 
the structure factor S{K,q = 2/L) ^ 0.1, while values 
for different q are close to 0, which happens to be the 
case for K sufficiently far from Kc. On the other hand, 
for DW, S{K,q = 1) = 0.25 and vanishes for other q 
values. Such a behavior is seen close to the resonance, 
K > 2.35. Thus indeed the two phases obtained close to 
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the resonance and far from it show the properties pre¬ 
dicted in [20]. Note that since the number of particles is 
strictly conserved in exact diagonalization, we cannot use 
some mean field order parameter to classify the phases 
observed. Still the identification based on the structure 
factor is unambiguous. 

The behavior is more complicated in the intermedi¬ 
ate interval of K values. The structure factor for both 
q = 2/L and q = 1 becomes small while intermediate 
q values {Aj..) become important. The situation 
seems somehow clearer close to the border of phase tran¬ 
sitions. Around K = 2.3 the peak in fidelity suscepti¬ 
bility coincides with the change in ground state struc¬ 
ture (as seen in S{K,q) plot) - instead of the fully sep¬ 
arated phases of composites and empty sites we observe 
splitting of the composites cluster into two (in small 
{AK ^ 0.001) interval directly above K = 2.3) and three 
clusters (which corresponds to the dominant S{K,3/7) 
value). Let us denote the pure CL phase as a string 
OOOOOOOCCCCCCC with C sites being filled by compos¬ 
ites. Respective many-cluster phases can be traced back 
to OOOOOCCOOCCCCC and OOOCCOOCCOOCCC configu¬ 
rations as verified by a careful examination of the ground 
state wave-function expansion in Fock space (possible due 
to the small size of our system). On the other hand while 
in the vicinity of Kc we observe a pure DW phase, close to 
K = 2.35 the inspection of the wave-function reveals an 
addition of defected component, with 2 sites breaking the 
DW symmetry. The relative importance of a single defect 
component changes smoothly from practically zero close 
to Kc (observe that above K = 2.35 all g-components 
of S{K,q) vanish except g = 1) to become significant 
below K = 2.35. The subsequent peaks of the fidelity 
susceptibility, y in Fig. 4 nicely coincide with different 
components of S{K,q) dominating the structure factor. 
That corresponds, as again confirmed by the inspection 
of the wave-function components, to successive defects 
of the partial DW leading to small clusters eventually 
merging as K moves further away from Kc. 

One may pose an important question whether the 
mixed phase observed is not really a finite size effect 
which will disappear in the thermodynamic limit and the 
mean field analysis [20] will be recovered in that limit. 
To provide an answer, we have evaluated the borders be¬ 
tween different phases for the longer chain with L = 18. 
Plotting the borders as a function of l/L and extrap¬ 
olating to the infinite chain one can see a clear indica¬ 
tion that the mixed phase should not be purely a finite 
size effect. Let us note that this behavior is reminiscent 
of striped phases observed in two-dimensional Falicov- 
Kimball model [39]. Importantly, considering the stan¬ 
dard optical lattices systems, the typical lattice size is 
about 50 sites thus the results obtained here are of a 
direct experimental relevance. 

Using results from diagonalizations, one can easily 
calculate the correlation function of composite bosons, 
(cgCj). When a system is in DW phase, the corre¬ 
lation function decays exponentially with increasing j 



FIG. 6. Borders between DW and MX phase as well as MX 
and CL phase as a function of the inverse system size. Observe 
that the mixed phase persists in the extrapolated L ^ oo 
limit. 
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FIG. 7. Correlation function (cjcj) of creation operators for 
boson composites in ground state calculated for L = 14 sites 
system for three values of parameter K — K/uj^ which are 
in the mixed {K — 2.32) and density wave {K — 2.36, 2.40) 
phases. The correlation function decays exponentially in the 
DW phase, for the MX and CL (not shown) phase the decay 
is much slower, presumably power-like. 


|(cQCj)| (X exp(—jf’//c) - compare Fig. 7. The correlation 
length Ic depends strongly on ^ - compare the correla¬ 
tion functions for K = 2.36 and K = 2.40. For other 
phases much slower decay, presumably power-like, is ob¬ 
served but no definite conclusions may be drawn due to 
small sizes considered. To that end, one should perform a 
numerical study of a much larger chain, e.g. using density 
matrix renormalization group (DMRG) which is beyond 
the scope of the present work. 
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FIG. 8. Different phases in the K - A plane. Observe the 
shrinking of the intermediate mixed phases region while the 
interesting density wave arrangement of composites region 
gets larger. The borders obtained for this plot were collected 
from diagonalizations of L = 10 system with 5 composites, 
results for L = 14 at selected points show that the picture is 
not affected significantly by the system size. 


B. Detuned case 


While the studies of [20] and the results presented 
above concentrated on the A = 0 case corresponding in 
DW phase to SSH Hamiltonian [ 12 ], it is interesting to 
see whether the full Rice-Mele model [13] for A 7 ^ 0 leads 
to similar conclusions. To that end, we have studied the 
phase diagram in the K — A plane as shown in Fig. 8 . 
Observe that while the border of the CL phase is not sen¬ 
sitive to A, the region of DW actually increases eating 
up the MX phase. Therefore, the Rice-Mele model seems 
to be realized quite easily with the present system. 

The most interesting physics of the Rice-Mele model 
comes from localized modes on defects on the borders be¬ 
tween topological and trivial phases [13, 40]. As discussed 
in [ 20 ], the present model allows for control of the num¬ 
ber of defects by changing slightly the filling of minority 
fermions, i.e. of composites. fFr the filling ric < 1/2 one 
creates holes in the DW, for ric > 1/2 we should have 
extra particles. Indeed, as visualized in Fig. 9 when we 
consider 6 composites in L = 13 sites, the DW phase (oc¬ 
curring for 7 composites and L = 14 sites) is replaced by 
a “single hole phase” (SHP). Due to the periodic bound¬ 
ary conditions and the translational invariance of the sys¬ 
tem, the ground state is a combination of states with a 
hole at different positions along the lattice as revealed 
by the eigenstate inspection in the Fock representation. 
The border between a SHP and mixed configurations is 
placed close to the value for the border of the DW phase 
in an ideal half filling of composites (taking into account 
finite size effects). For K far from Kc^ we observe a sharp 
phase transition to clustered phase with holes and defects 
separated. 


FIG. 9. Structure factor S{K,q) for L = 13 and N = 6 
composites. Such a situation results in a defect (hole) present 
in the system. A single hole leads to a structure factor being 
dominant for q = 12/13, the presence of a hole results in non 
vanishing values of S for other g’s. As before, the GL phase 
is characterized by the structure factor being most prominent 
at q = 2/L. 


IV. CONCLUSIONS 

Using exact diagonalization on small systems, we have 
addressed the problem of resonantly shaken optical lat¬ 
tices in which an unevenly populated mixture of two 
species of fermions is held. We have verified the basic 
model studied in [ 20 ] where, neglecting minority fermion 
tunnelings, a density wave arrangements of composites 
were found in the situation when the shaking ampli¬ 
tude was tuned in a way enabling switching off all of 
the intra-band tunnelings. Then, the excess majority 
fermions move in an emergent lattice (formed by com¬ 
posites) with direction dependent tunnelings realizing the 
topological Rice-Mele model. In the simplified approach 
[ 20 ], it was found that apart from the density wave (for 
switched off intra-band tunnelings) the composites and 
empty sites may separate forming two clusters - when 
the intra-band tunnelings are important. That has also 
been confirmed by the present calculation. In addition 
to these two phases, the middle region separating these 
ideal case is revealed by an exact diagonalization. In this 
mixed-phase region, the ground state contains superpo¬ 
sition of many different composite arrangements. This 
phase may show quasi-long range order which is absent 
in the density wave phase. 

We have also shown that the density wave phase in the 
vicinity of shaking parameters combination switching off 
intra-band tunnelings {K/uj ^ 2.405 - the zero of 0-order 
Bessel function) persists even when the shaking frequency 
is not adapted precisely to the s-p orbital resonance con¬ 
dition - thus it is quite robust. We have explicitly shown 
that the deviations from the ideal half filling of the mi¬ 
nority fermions (and thus the composites) leads directly 
to defects (holes or extra particle) that, if occurring on 
the edge of the topologically nontrivial phase lead to lo- 
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